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INTRODUCTION 


One  of  the  most  serious  limitations  of  single  scattering  code  AGAUS,  as 
developed  under  previous  contracts,  was  the  Inability  of  the  Internal  Mle 
routine  to  provide  reliable  results  for  large  diameter  particles  and 
situations  Involving  relatively  large  absorption.  These  limitations  have  been 
reduced  substantially  through  the  replacement  of  the  older  forward  recursion 
Mle  routine  by  one  using  partial  backward  recursion  and  a  method  of  continued 
fractions. 

The  new  routine  has  been  made  operational  on  the  HP  computer  system  at  the 
US  Army  Atmospheric  Sciences  Laboratory,  White  Sands  Missile  Range,  NM,  and 
has  been  coded  to  serve  as  a  direct  substitute  for  the  earlier  routine.  This 
report  has  been  written  to  serve  primarily  as  documentation  of  the  new 
routine.  For  the  sake  of  completeness  a  brief  review  of  the  Mle  theory  has 
been  Included  as  well  as  symbolic  definitions,  flow  charts  for  the  new  routine 
and  a  discussion  of  Its  reliability. 

It  should  be  noted  that  the  new  routine  has  been  coded  to  have  the  same 
name  (MIEGX)  as  the  older  Mle  routine.  This  was  done  for  ease  of  replacement 
by  users  who  do  not  wish  to  revise  other  portions  of  AGAUS  to  reflect  a  change 
of  name.  The  new  routine  can  be  easily  distinguished  from  the  older  one  by 
Its  use  of  complex  arithmetic. 


5  paos  bunk-not 


MIE  THEORY* 

Mie  theory  predicts  the  scattering  by  and  the  absorption  in  an 
isolated,  discrete,  homogeneous,  isotropic  sphere  of  diameter  D  with 
a  known  complex  refractive  index  n  =  m-ik  relative  to  the  surrounding 
medium  and  illuminated  by  monochromatic  radiant  energy  with  wavelength 
X  in  the  surrounding  medium.  The  theory  is  given  in  detail  in  standard 
texts  and  need  not  be  repeated  here.  Instead,  those  elements  of  theory 
needed  for  an  understanding  of  the  numerical  algorithms  used  are 
included. 

All  scattering  properties  of  spheres  are  computed  from  m  and  k, 

and  through  the  use  of  the  induced  electric  and  magnetic  multipole 

moments  of  the  sphere  a  and  b  ,  respectively.  The  moments  are  given  by^ 

n  n 


'f^(na)  H'^(a)-n  'i'^(na)  r  (a) 

^n  "  i"(na)  C  (a)-n  (n  )  C’ (a)'  ’ 
n  n  n  n 

and 

n  'f'  (na)  Y  (a)-Y  (no)  Y'  (a) 

,  _  n _ n _ n _ n _ 

n  n  Y'(na)  5  (a)-Y  (no)  ^'(a) 
n  n  n  n 


(1) 

U) 


This  section  is  partially  based  on  material  taken  from  ECOM  report 
ECOM-5558  by  R.B.  Gomez,  C.  Petracca,  C.  Querfeld  and  G.B.  Hoidale, 
March  1975,  and  the  final  report  on  Contract  DAAD07-78-C-0063  by  Miller 
et  al. ,  December  1978. 


Note  that  n  is  used  as  a  subscript,  an  integer  index  and  a  complex 
index  of  refraction  when  it  is  not  a  subscript. 
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The  prime  denotes  differentiation  with  respect  to  the  argument.  The 

Y  (z)  and  5  (z)  functions  are  Ricatti-Bessel  functions  of  the  first 
n  ^ 

arid  third  kind,  respectively,  and  are  related  to  spherical  Bessel 

functions  j  (z)  and  n  (z)  by 
n  n 


.  z  J^(z), 

and 

5^(2)  =  zj^(z)-i  zn^Cz)  =  'l'^(z)  +  ix^(z), 

where 


(3) 


(4) 


(5) 


and 


n  (z) 


^2z^  \+l/2^^^' 


(6) 


The  function  half  integral  order  Bessel  function;  the 

function  ^^+1/2^^^  half  integral  order  Neuman  function. 

The  extinction  cross  section  is  computed  from 


<VV-  ('> 

n=l 

and  the  scattering  cross  section  from 

n=l 

The  various  cross  sections  are  the  basic  quantities  used  in  scattering 
problems,  but  they  are  not  the  quantities  usually  computed  directly  from 
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Mle  algorithms.  Instead,  it  is  more  convenient  to  compute  dimensionless 


efficiency  factors  Q  and  Q  ,  which  depend  on  n,  k,  and  a,  and  which 

6Xt  SC3 

are  multiplied  by  the  geometrical  sphere  cross  section  to  obtain  the  true 

cross  section  =  irr^Q . .  Thus, 

1  1 

CO 

Q  ,  =  -T  y  (2n+l)  Re(a  +b  ),  (9) 

ext  ^  '  n  n 

n=l 

and 

■isca  -  -if  ^ 

n=l 

Although  the  cross  sections  account  for  the  energy  removed  from  the 
forward  beam,  they  do  not  give  any  information  about  where  the  scattered 
photons  go.  This  information  is  contained  in  scattering  amplitudes  and 
intensity  factors  which  relate  the  flux  density  scattered  through  an 
angle  0  relative  to  the  incident  flux  density.  There  are  two  amplitudes, 
S^(0)  and  82(6) »  and  intensity  factors  lj^(0)  and  12(8),  which  correspond 
to  light  respectively  polarized  perpendicular  and  parallel  to  the  plane 
of  scattering  defined  by  the  direction  of  incidence  and  the  direction  of 
scattering. 

The  intensity  factors  are  related  to  the  scattering  amplitudes  by 


ij^(e)  =  |Sj^(0)|2  , 

(11) 

12(e)  -  |s2(e)p  , 

(12) 

1^(0)  =  ReCS^-S*},  and 

(13) 

i^(0)  =  -Imag{Sj^*S2}  ' 

(IM 
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The  amplitudes  come  from  the  multipole  moments  through 


=  I  -^T^T  [a  ir  (6)  +  b  T  (e)  ] , 
1  n(n+l)  n  n  n  n 

n=i 


5,(9)  =  I  TT  (6)  +  a  T  (0)], 

I  n(n+l)  n  n  n  n 

n=i 


and  angular  factors  tt  (.0)  and  t  (0J  defined  in  terms  of  associated 

n  n 

Legendre  functions: 


IT  (0)  =  P ^  (cos0) /sin0 , 
n  n 


dP^ (cos6) 


t^(0)  = 


The  amplitudes  have  relative  phase  6  =  argS^  -  argS^. 
Alternative  expressions  frequently  used  are 


dP  (cosO) 
n _ 

d (cosO) 


dir  (0) 

ir  (6)  =  cose-rr  (0)  -  sin^e  .  —— — ^  , 
n  n  dlcose) 


where 


P  (cos0) 
n 


1  d 

2^n!  dcos^S 


-(cos^6-l)”. 
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These  functions  satisfy  the  following  recurrence  relations: 


71  (6)  •=  cose  W  ,  (e)  -  77  (6), 

n  (n-1)  n-1  n-1  n-2 


(22) 


and 


T  (6)  =  cosBItt  (6)-7I  „(6)  ]-(2n-l)sin^6iT  ,(6)+77  ~(e).  (23) 

n  n  n-2  n-1  n-2 


The  scattering  cross  section  measures  the  ability  of  a  particle  to 
scatter  light,  and  it  is  to  be  expected  that  C  is  obtained  from  an 

SC  d 

integral  over  the  scattering  intensity  factors.  Equation  (8)  follows  from 


C 

sea 


1 

(i^(e)  +  12(6))  dcose. 

-1 


(24) 


Although  the  intensity  factors  themselves  may  be  used  in  scattering 
calculations,  they  are  primarily  suited  for  computing  flux  densities, 
and  it  is  frequently  more  convenient  to  measure  and  compute  scattered  light 
in  tenns  of  radiances.  Radiances  do  not  have  the  1/r^  dependence,  and  it 
is  therefore  unnecessary  to  know  the  distance  from  the  scatterer  to  the 
detector  if  the  detector  field  ot  view  is  small  and  is  filled  by  the 
scattering  cloud.  The  phase  function  pi0)  gives  a  radiance  I  scattered 
into  the  9  direction  in  terms  of  the  radiance  incident  on  the  particle. 

The  phase  function  is  dimensionless  and  is  defined  here  as 


p(e)  =  2-^— 

ext 

The  normalized  phase  function  p(6)dn/477  gives  the  probability  of  a 
photon  being  scattered  through  an  angle  6  into  an  element  of  solid 
angle  dfi  *  difdcosB.  The  integral  of  the  normalized  phase  function  is 
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(33) 


SUBROUTINE  MIEGX 

Subroutine  MIEGX  computes  various  efficiency  factors,  and  intensity 
factors  i^,  each  complex  refractive  index  m  and  size 
parameter  a  (or  x) .  The  Ricatti-Bessel  functions  and  their  derivatives 
in  Eqs.  (1)  and  (2)  are  computed  by  the  forward  recursion  method. 


£  (z)  =  £  ,(z)  -  £  Az),  (37) 

n  2  n-i  n-2 

where 

£^^(2;)  =  ^  complex  quantity. 

The  initial  values  used  in  forward  recursion  are 


12 


(2)  =  sin  2, 
o 

,  ,  ^  sin  2 

(z)  - - cos  2 , 

I  2 


X^(2)  =  COS  2,  and 


cos  2 

X  =  -  +  sin  2. 


The  angular  functions  tt  and  t  are  also  computed  by  forward  recursions 

n  n 

from  Eqs.  (22)  and  (23).  The  initial  values  used  are  it  (6)  =  0,  it,  (6)  =  1, 

o  1 

T  (6)  =  0,  and  t  (6)  =  cos6. 
o  1 

For  computational  purposes  it  is  more  convenient  to  write  Eqs.  (1)  and  (2) 
+ 


in  the  following  forci 


A  (n  a) 
n  1  1  n 

Re(C  (a)]  -  ReIC 

,(a)] 

n^  a 

n 

n-1 

n 


A  (n.a) 
n  1 


C  (a)  -  E  ,(c.) 
n  n-1 


(38) 


[n  A  (n  a)  +  ^]  Re[C  (a)]  -  Re[C  ,  (a) ] 
b  =  — UL-i - 2 - n - n::! -  ^ 

"  In  A  (n  o)  +  ^]  C  (a)  -  E  (a) 

1  n  1  an  n-1 


where  A  (n.a)  = 
n  1 


J  ,  (n  u) 
_  H  +  v-i  1 

z  J  (n.a) 

V  i 


(^0) 


and  V  =  n  +  -J.  The  symbol  n^  is  used  here  for  the  complex  index  of 

refraction  to  distinguish  it  from  subscript  n. 

Methods  utilizing  either  forward  or  backward  recursion  on  the  ratio 

have  been  applied  to  the  problem  and  both  methods  have  their  own  unique 

problems.  Forward  recursion  methods  are  limited  by  the  number  of 

+"Development  of  Programs  for  Computing  Characteristics  of  Ultraviolet 
Radiation,"  IBM  Corp.,  1972,  Technical  Report  for  Contract  No.  NASS-21680 
(NASA) 
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significant  digits  in  the  computer.  Backward  recursion  methods  must, 

of  course,  calculate  for  some  N  larger  than  the  n  required  for 

convergence  of  the  Mle  sums.  Both  methods  usually  require  double 

precision  arithmetic  and  have  been  shown  to  fail  for  cases  Involving 

large  a  and/or  a  large  Imaginary  part  of  the  refractive  index. 

This  subroutine  uses  the  method  of  continued  fractions^  to 

calculate  the  ratio  independent  of  any  previous  value.  The  ratio 

is  correct  to  the  accuracy  of  the  machine.  The  values  of  A  for  n<N 

n 

are  then  calculated  using  the  backward  recursion  formula 


2n-l 

z 


.i  Cz) 


(^1) 


If  convergence  of  the  I'fie  summations  requires  n>N,  then  A2j^  is 

calculated,  again  independent  of  previous  A  *s,  and  backward 

n 

recursion  is  used  to  calculate  A  for  N+l<n<2N. 

n 


This  description  is  based,  in  part,  on  technical  report  ECOM  5509  (1973) 
AD767223  by  W.J.  Lentz.  An  article  covering  the  same  material  can  be 
found  in  Applied  Optics  No.  3,  b68  (1976). 
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A  continued  fraction  may  be  written  as 


f  -  ao  + 


c_^  +  b2 


^2  ^3 


c  ^  •  •  • 


f  =  +  ^  ^  ^ 

^0  c  1"  c  c  1“ 

1  2  3 


The  n  approximate  convergent  to  the  continued  fraction  representation 
of  f  is  written  as 


b  b 

,  _  ,  n 

n  ^0  c  ^  c 

1  n 


Continued  fractions  may  be  generated  from  a  three-term  recursion  relation 
* 

in  a  simple  way: 


2(v+l)  _  v+2 

X  J  .  . 


_  1  1 
z  ~  2(v+l)  _  2(vf2) 


=  2vz 


2(v+i)z-I  -  2(^^^-2)z“l  - 


For  simplicity  the  argument  z  may  be  suppressed. 
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This  form  is  one  of  a  simple  continued  fraction  which  may  be  defined  as 


a  +  —  —  — 

1  32  +  +  . . 


where  a^^  0  and  the  a's  may  be  negative.  Equation  (47)  may  be  written 

in  the  more  convenient  notation  of 


f (x)  —  [a^^  ^2*  ^3*  • * • 3 • 


The  n  convergent  is  written  in  like  manner: 


fn(x)  =  [a^,  32,  a^,  ...  ,  a^] 


Lentz  has  shown  that 


[a^]  ...  [a  .< 

f  (x)  ®  7 - ; - f - 

n  ^^2*^  •••  .2. *  ** 


•  *  •••  » 

.  ,  a2^^®n’  ***  ’  ^2^ 


The  advantage  of  this  method  is  that  the  calculation  begins  with  the  first 
term  (rather  than  the  end)  of  the  fraction  and  is  terminated  when  it  is 
determined  that 


[aj^,  • '  •  »  a^^] 


-  1  <  e  . 


The  value  of  e  depends  on  the  accuracy  needed  or  on  the  number  of 

significant  digits  available  on  the  computer  used. 

The  ratio  J  ,/J  can  be  written  in  a  similar  form  with 
v-1  V 


a  -  (-1)"  2(v+n+l)z”^. 

n 


Subroutine  MIEGX,  as  furnished,  terminates  the  calculation  of  A 

n 

with  e  =  10“^.  This  could  be  decreased  to  increase  speed  at  the  loss 

of  some  accuracy.  If  the  routine  was  converted  to  double  precision 

it  would  be  possible  to  set  e  to  a  much  smaller  value. 

MIEGX  then  computes  the  preceding  values  of  A  by  baclcward  recursion 

n 

and  stores  them  to  be  used  in  the  Mie  summations  of 

Re[S^(e)],  Im[S^(e)],  Re[S2(e)]  and  Im[S2(e)]. 

The  sum  is  terminated  when 

and  when  the  fractional  change  in  the  radar  efficiency  is  also  less  than 
c ,  1 . e . ,  when 

^rad  n  ^rad  n-1 

_  ^  c  • 

^rad 

This  is  more  stringent  than  the  first  test  alone  and  is  a  test  on  the 
phase  functions  as  well. 

MIEGX  returns  the  following  quantities  as  required  by  AGAUSX: 

Q  Q  .  Q  P(J).  PFNZRO,  OISTAR,  and  02STAR.  The  P(J)  and  PFNZRO 
6XC  scd  irso 

as  returned  by  MIEGX  are  average  intensities  and  must  be  further  normalized 
to  become  the  actual  phase  functions  (Eq.  25) .  OISTAR  and  02STAR  are  the 
Legendre  expansion  coefficients  and  co^* 


VALIDATION  OF  THE  CODE 


Comparison  runs  ot  AGAUS  using  MIEGX  with  continued  fractions  and 
a  previous  version  of  MIEGX  based  on  forward  recursion  were  done  in 
the  region  that  the  forward  recursion  routine  is  considered  valid  (see 
final  report  on  Contract  DAAD07-78-C-0063,  December  1978).  The 
comparison  process  was  done  mainly  to  verify  that  coding  errors  had  been 
eliminated  and  that  interfacing  of  the  two  was  accomplished  properly. 

When  one  begins  to  test  the  limits  of  a  code  it  is  important  to 

keep  in  mind  the  predictions  of  Mie  theory.  In  the  limit  as  a 

is  expected  to  converge  to  2.0  for  constant  m  and  k.  There  may  be 

fluctuations  or  ripples  but  these  will  be  small  for  large  a.  The  absorption 

efficiency  factor,  is  dependent  on  the  imaginary  part  of  the  index 

of  refraction.  will  be  small  for  the  small  imaginary  part  and  approaches 

1  for  the  totally  absorbing  sphere.  Qg^.^  is  bhe  difference  between 

Q  and  Q  .  and  therefore  is  bounded  above  by  Q  .  for  small  imaginary 
ext  ^abs  ext 

part  and  equal  to  for  an  all  real  index  of  refraction.  For  complex 

indices  will  approach  a  value  between  Qgxt~^abs  ^ext’  i'®’» 

between  1  and  2  as  a  becomes  large  but  at  no  time  may  exceed 

The  radar  efficiency  factor  is  expected  to  approach  a  limit  between  0  and 
1  as  a  increases.  The  larger  the  value  of  k  the  smaller  the  asymptotic 
value  of  i®  expected  to  be.  Some  oscillation  is  expected  and  it 

is  possible  that  resonances  can  produce  peak  values  much  larger  than  the 
asymptotic  value  and  minima  near  zero.  Resonances  may  be  more  pronounced 
for  k  approximately  0.  Values  of  larger  than  10  would  probably  be 

in  error. 
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Table  1  shows  the  results  of  comparison  calculations  using 

MIEGX  and  DAMIE  (a  Mle  routine  written  by  J.  V.  Dave).  The  values  of 

m,  k,  and  a  presented  were  chosen  to  show  the  agreement  between  the  two 

codes  as  well  as  to  point  out  the  extended  usefulness  of  MIEGX.  0  . 

^radar 

was  included  in  the  comparisons  because  it  appears  to  be  very  sensitive 
to  computational  errors  since  it  depends  on  the  difference  between  a 

n 

and  b  .  Q  .  is  also  an  indicator  of  the  validity  of  the  phase  functions 
n  rsuHr 

being  Itself  p(180*). 

For  n  =  1.2  -  O.Oi  and  a  =  10  to  400, there  is  exact  agreement  between 

MIEGX  and  DAMIE  on  the  value  of  Q  ^  and  Q  .  There  is  some  variation 

ext  ^sca 

between  the  two  for  Most  interesting,  though,  is  the  large 

result  for  a  =  50  obtained  with  both  routines.  This  might  be  explained 

by  resonance  as  mentioned  previously. 

The  results  obtained  for  n  =  1.2  -  0.6i  show  the  failure  of  the 

DAMTE  routine  for  large  a.  As  discussed  above  Q  must  be  less  than 

sea 

A  value  of  7.5  (a  =  200)  is  definite  proof  of  failure.  For  a  =  100 

and  larger,  deviated  dramatically  from  the  value  obtained 

with  MIEGX.  The  presence  of  a  slow  downward  trend  in  the  MIEGX  values 

of  Q  ^  and  Q  as  well  the  near  constant  value  of  Q  ,  indicate  that 

ext  sea  radar 

the  routine  is  probably  still  valid  for  a  of  400.  The  results  for 
n  »  1.2  -  1.21  are  presented  to  again  display  the  extended  usefulness  of 
MIEGX  with  continued  fractions. 

Based  on  these  results,  as  well  as  other  test  results  not  Included 
in  this  report,  MIEGX  is  likely  to  be  reliable  for  a’s  of  at  least  400 
with  an  index  of  refraction  as  large  as  5-51.  For  the  wide  variety  of  test 
cases  using  MIEGx  there  has  never  been  an  obvious  error  in  the  calculation 
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If  the  user  has  an  application  that  requires  unusual  sizes  and/or 
index  of  refraction  combinations  It  is  suggested  that  test  runs  be 
done  that  encompass  the  region  of  interest.  Comparing  the  results  to 
theory  should  provide  an  Idea  of  the  accuracy  In  that  region. 
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Table  1.  Ce.barl.oe  ot  KIEOC  (vltb  Ceetlaeed  Fractleaa)  aiai  MMb 


List  of  Symbols  used  In  MltCX 


SYMBOL 

A(J) 

ACAPN 
AXPHA.X 
ALPHAD 
C(J) ,CJ 

CAV.CAYD 

CAYE 

ELTRHX(1,J) 

ELTRMX(2,J) 

ELTRMXO.J) 

ELTRMX(4,J) 

EM. EMD 

EN. ENLl 

FNA, FNAP,FNAPP 

FNB, FNBP,FNBPP 
N 

NDELTA 

NDIM 

NMX 

NMIN 


_ Explanation  or  Definition _ 

A  (Eq.  40)  the  array 
n 

A  for  the  current  n  of  the  Mle  sum 
n 

Mie  size  parameter,  a  =  X  =  2vr/\ 

Double  precision  ALPHA 

The  array  of  cosines  of  the  scattering  angles;  there  are 
'it'  elements  in  the  array,  the  Jth  element,  respectively. 

The  ratio  of  the  imaginary  part  to  the  real  part  of 
adjusted  refractive  index,  double  precision  CAY,  respectively. 

The  true  Imaginary  part  of  adjusted  refractive  index 

Within  n  loop:  =  Re[S^(0)];  after  n  loop:  =  I2  = 

S2(e)-S*(e) 

Within  n  loop:  =  Im[S^(6)];  after  n  loop:  =  1^  = 

S^(6)-S*(e) 

Within  n  loop:  =  RefS2(6)];  after  n  loop:  =  ® 

Re[Sj^(e)-S2(e)] 

Within  n  loop:  =  Im[S2(6));  after  n  loop:  =  ® 

Im[S^(0)S2(e)] 

The  real  part  of  adjusted  refractive  index,  double 
precision  EM,  respectively. 

Floating  point  representation  of  N,  N-1 ,  respectively. 

=  a  ,  a  ,,  a  _,  respectively, 
n  n-1  n-2 

=  b  ,  b  ,,  b  _,  respectively, 
n  n— 1  n-z 

Index  in  Mie  sum 

The  smaller  of  NDIM  and  NMX,  later  used  as  increment 
of  N  for  calculation  of  A 

n 

The  dimension  of  A  array 

“  X*(nH-k)+9  an  approximation  to  the  maximum  N  needed 
-  NMX+1 -NDELTA 
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SYMBOL 


Explanation  or  Definition 


OISTAR.OISTRD 

U2STAR,02STRD 

P(J) 

PFNZRO 

PICl  .J)  .PI(2.J) 
PI(3,J) 

PI1J,PI2J,PI3J 

QEXT 

QRD 

QRT 

QRTL2 

QRTR 

QSCAT 

QSD 

QTD 

RE 

RRF 

RRFX 

RX 

SGR 

SGS 

SGT 

SUMRR 


=  first  order  coefficient  for  Legendre  expansion  of 

the  average  intensity  P(J),  douole  precision  OISTAR  respectively 

-* 

=  “2 

=  (i^+i2)/2,  tne  average  intensity  at  angles,  arc  cos(C(J)). 

'IT'  elements  in  array 

=  the  average  intensity  at  0“ 

=  IT  -(e),  tr  .(6),  Tt  (0)  (Eq.  22) 
n-2  n-1  n-j 

The  Jth  element  of  PI(1,J),  PI(2,J)  and  PI(3,J), 
respectively 

Unnormalized  efficiency  factors 

Double  precision  backscatter  efficiency  factor 

=  SUMRR^  +  SUMRI^  present  value  of  unnormalized 
backscatter  cross  section 

Previous  value  of  QRT 

=  I  QRT  -  QRTLl|/QRT,  ratio  of  change  in  QRT  to  present 
value  -  used  as  part  of  exit  criterion 

The  unnormalized  scattering  cross  section 

Double  precision  scattering  efficiency  factor 

Double  precision  extinction  efficiency  factor 

=  EM  -  iCAYE,  complex  refractive  index 

=  1/RF 

=  1/(X*RF) 

=  1/X 

The  backscattering  (radar)  efficiency  factor 

The  scattering  efficiency  factor 

The  extinction  efficiency  factor 
N 

»  y  (-l)"(2n+l)  Re [a  -b  ] 

^  '  n  n 

n=l 
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TE(1) ,TE(2J 
TF(1) ,TF(2) 
TG(1),TG(2) 
TCI 

TC2 

TOL 

V 

vmi 

WFN(l) 

WFN(2) 

X,  ALPHA 

Y 

ZAN.ZANP 

ZDEN.ZNUM 


Real  and  Imaginary  part  of  FNBP,  respectively 

Real  and  imaginary  parts  of  FNAPP,  respectively 

Real  and  Imaginary  parts  of  FNBPP,  respectively 

=  (A  /n.)  +  Cn/x) 
n  1 

=  (A  ‘n  )  +  (n/x) 
n  1 

«  l.E-06  exit  tolerance  for  Mie  sum 
•=  NMX  +  3/2 


-  X*RF 


Temporary  complex  valued  variables  in  continued 
fractions  calculation  of  A 

n 
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MIEGX  -  Simplified  Flowchart 


NOTES  ON  FLOKCiLVRT  RATIO  OF  CONSECUTIVE  BESSEL  FUNCTION 


Wliat  seems  to  be  a  confusing  use  of  variables  within  the 

calculation  loop;  namely  NMX,  NMIN,  NDELTA,  JJ  and  J,  is  necessary  because 

the  wav  the  A  array  is  used.  The  first  time  through,  the  ratio  J  ,/J 

and  .Ajj  is  calculated  for  the  largest  value  of  N  that  might  be  needed  or 

is  allowed  by  array  size.  The  ratio  is  calculated  using  the  continued 

fractions  routine  and  the  preceding  A^'s  are  calculated  using  backward 

recursion.  If  the  Mie  sum  does  not  converge  before  it  uses  all  the 

calculated  A  then  it  becomes  necessary  to  calculate  the  next  A  values, 
n  n 

If  so  the  ratio  and  A^,  are  recalculated  for  a  new  largest  N  (this  time 
twice  the  original  N)  and  the  previous  N  values  of  A^  replace  the  former 
array  elements,  l.e.,  the  A  array  now  contains  •••  ^2N’  this 

is  done  the  Mie  sum  continues.  This  calculation  and  replacement  of  the 
A  array  continues  until  the  Mie  sum  converges. 
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